import numpy.linalg as la
import numpy as np

ndim = np.array([2,3,8,11,14])

Let's perform linear solves for matrices with increasing size "n", for a problem in which we know what the solution would be.

for nd in ndim:
    ## This is the vector 'x' that we want to obtain (the exact one)
    x = np.ones(nd)
    ## Create a matrix with random values between 0 and 1
    A = np.random.rand(nd,nd)
    ## We compute the matrix-vector multiplication 
    ## to find the right-hand side b
    b = A @ x
    ## We now use the linear algebra pack to compute Ax = b and solve for x
    x_solve = la.solve(A,b)
    ## What do we expect? 
    print("------ N =", nd, "----------")
    error = x_solve-x
    print("Norm of error = ", la.norm(error,2)) 

------ N = 2 ----------
Norm of error =  0.0
------ N = 3 ----------
Norm of error =  1.1957467920563633e-15
------ N = 8 ----------
Norm of error =  2.1339110374783316e-14
------ N = 11 ----------
Norm of error =  3.7696281795943905e-14
------ N = 14 ----------
Norm of error =  2.0764138859057344e-14

print(x_solve)

[1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]

Now we will perform the same computation, but for a special matrix, known as the Hilbert matrix

def Hilbert(n):
    
    H = np.zeros((n, n))    
    for i in range(n):        
        for j in range(n):        
            H[i,j] = 1.0/(j+i+1)    
    return H

for nd in ndim:
    ## This is the vector 'x' that we want to obtain (the exact one)
    x = np.ones(nd)
    ## Create the Hilbert matrix
    A = Hilbert(nd)
    ## We compute the matrix-vector multiplication 
    ## to find the right-hand side b
    b = A @ x
    
    ## We now use the linear algebra pack to compute Ax = b and solve for x
    x_solve = la.solve(A,b)
    ## What do we expect? 
    print("------ N =", nd, "----------")
    error = x_solve-x
    print("Norm of error = ", la.norm(error,2)) 

------ N = 2 ----------
Norm of error =  8.005932084973442e-16
------ N = 3 ----------
Norm of error =  1.389554002205336e-14
------ N = 8 ----------
Norm of error =  5.035065902200173e-07
------ N = 11 ----------
Norm of error =  0.019591809632407014
------ N = 14 ----------
Norm of error =  6.999747915358066

print(x_solve)

[ 0.99999985  1.00002076  0.99927926  1.01069926  0.91619132  1.38032247
 -0.02036746  2.47417863  0.56641165 -1.35352445  5.55922945 -2.98137159
  2.77700163  0.67192909]

What went wrong?

Condition number

The solution to this linear system is extremely sensitive to small changes in the matrix entries and the right-hand side entries. What is the condition number of the Hilbert matrix?

for nd in ndim:
    ## This is the vector 'x' that we want to obtain (the exact one)
    x = np.ones(nd)
    ## Create the Hilbert matrix
    A = Hilbert(nd)
    ## We compute the matrix-vector multiplication 
    ## to find the right-hand side b
    b = A @ x
    ## We now use the linear algebra pack to compute Ax = b and solve for x
    x_solve = la.solve(A,b)
    ## What do we expect? 
    print("------ N =", nd, "----------")
    error = x_solve-x
    print("Norm of error = ", la.norm(error,2)) 
    print("Condition number = ", la.cond(A,2))

------ N = 2 ----------
Norm of error =  8.005932084973442e-16
Condition number =  19.281470067903975
------ N = 3 ----------
Norm of error =  1.389554002205336e-14
Condition number =  524.0567775860629
------ N = 8 ----------
Norm of error =  5.035065902200173e-07
Condition number =  15257574847.190956
------ N = 11 ----------
Norm of error =  0.019591809632407014
Condition number =  521712392909565.75
------ N = 14 ----------
Norm of error =  6.999747915358066
Condition number =  6.980917584970002e+17

Residual

for nd in ndim:
    ## This is the vector 'x' that we want to obtain (the exact one)
    x = np.ones(nd)
    ## Create the Hilbert matrix
    A = Hilbert(nd)
    ## We compute the matrix-vector multiplication 
    ## to find the right-hand side b
    b = A @ x
    ## We now use the linear algebra pack to compute Ax = b and solve for x
    x_solve = la.solve(A,b)
    ## What do we expect? 
    print("------ N =", nd, "----------")
    error = x_solve-x
    residual = A@x_solve - b
    print("Error norm = ", la.norm(error,2)) 
    print("Residual norm = ", la.norm(residual,2)) 
    print("Condition number = ", la.cond(A,2))

x_solve

Rule of thumb

for nd in ndim:
    ## This is the vector 'x' that we want to obtain (the exact one)
    x = np.ones(nd)
    ## Create the Hilbert matrix
    A = Hilbert(nd)
    ## We compute the matrix-vector multiplication 
    ## to find the right-hand side b
    b = A @ x
    ## We now use the linear algebra pack to compute Ax = b and solve for x
    x_solve = la.solve(A,b)
    ## What do we expect? 
    print("------ N =", nd, "----------")
    error = x_solve-x
    residual = A@x_solve - b
    print("Error norm = ", la.norm(error,2)) 
    print("|dx| < ", la.norm(x)*la.cond(A,2)*10**(-16))
    print("Condition number = ", la.cond(A,2))